function [xdot,thrust] = rocketODEThermo(t,x,data)
% function to return derivatives for water bottle rocket state equations
% inputs are t = time (scalar), x = state vector (see below), and data
% (rocket adn atmospheric conditions)
%
% necessary components in data structure are:
% data.tt, data.thrust -> time and thrust vectors from static test data,
% interpolated, if no thrust data (solely glide, such as necessary for ISP
% approximations, initialize with tt = [0,1000], thrust = [0,0]
%
% data.wind -> 3-component vector of wind velocity
% data.initH -> normalized vector of initial heading when within launch pad
% data.bottleMass -> mass of empty rocket
% data.airDens -> density of air (atmosphere)
% data.bottleXArea -> rocket cross-sectional area for drag computation
% data.CD -> rocket Coefficient of drag
% data.gavec -> 3-component vector of gravitational acceleration 
% data.rhoW -> density of water or 'rocket fuel'....
% data.nozzleArea -> cross sectional area of the nozzle/bottle neck/exit
% data.bottleVol -> total volume of bottle
% data.initPressure -> fill pressure of bottle (Pa)
% data.coeffDischarge -> discharge coefficient of nozzle
% data.gamma -> specific heat ratio
% data.atmPress -> atmospheric pressure
% data.RealGasConst -> real gas constant for air
% data.ambientTemp -> ambient tempt of water/gas/everything


% x vec is [x,y,z,vx,vy,vz,mass water, mass air]


% work out thermodynamics
wVol = x(7)/data.rhoW;
airVol = data.bottleVol - wVol;
v0Air = data.bottleVol - (1/data.rhoW)*data.initWMass;
m0Air = data.initPressure*v0Air/(data.RealGasConst*data.ambientTemp);

if wVol > eps
    

    pair = data.initPressure*(v0Air/airVol)^(data.gamma);

    mdotW = -data.coeffDischarge*data.nozzleArea*...
        sqrt(2*data.rhoW*(pair - data.atmPress));
    
    mdotA = 0;
    
    thrust = 2*data.coeffDischarge*(pair - data.atmPress)*data.nozzleArea;
else
    % 'end' values are from the end of the previous stage, so the air
    % volume is that of the total volume
    pend = data.initPressure*(v0Air/data.bottleVol)^data.gamma;
    
    
    pair = pend*(x(8)/m0Air)^data.gamma;
    
    rho = x(8)/data.bottleVol;
    
    Temp = pair/(rho*data.RealGasConst);
    
    pstar = pair*(2/(data.gamma+1))^(data.gamma*(data.gamma-1));
    
    if pair - data.atmPress < eps
        ve = 0;
        rhoe = 0;
        pe = data.atmPress;
        
    elseif  pstar > data.atmPress %choked flow (Me = 1)
        ve = sqrt(data.gamma*data.RealGasConst*data.ambientTemp);
        Te = (2/(data.gamma+1))*Temp;
        pe = pstar;
        rhoe = pe/(data.RealGasConst*Te);
        
    else
        pe = pair;
        Me = findMachNumber(pair/data.atmPress,data.gamma);
        
        Te = Temp*(1+(data.gamma-1)/2*(Me^2));
        rhoe = pair/(data.RealGasConst*Te);
        
        ve = Me*sqrt(data.RealGasConst*data.gamma*Temp);
    end
    
    thrust = data.coeffDischarge*rhoe*data.nozzleArea*ve^2 + (pe-data.atmPress)*data.nozzleArea;
    mdotA = -data.coeffDischarge*rhoe*data.nozzleArea*ve;
    mdotW = 0;
end


 
% magnitude of position vector
posMag = norm(x(1:3));

% calc booleen as to whether rocket is still within pad (within 1m)
withinPad = posMag < 1;

% calculate relative velocity (take into account wind)
vrel = x(4:6) - data.wind;

%magnitude of wind relative velocity
vrelMag = norm(vrel);

%calculate heading, uses initial heading if within pad and the relative
%current velocity if not within pad
h = data.initH*withinPad + ~withinPad*vrel/(vrelMag + 1e-13);

totalMass = data.bottleMass + x(7) + x(8); %total mass of rocket is the bottle mass in addition to the water mass

%calculate drag magnitude
Drag = 0.5*data.airDens*(vrelMag^2)*data.bottleXArea*data.CD;

%calculate gravity force
gvec = data.gavec*totalMass;

%% lets output stuff
xdot(1:3) = x(4:6);

xdot(4:6) = (thrust*h - Drag*h + gvec)/totalMass;

xdot(7) = mdotW;

xdot(8) = mdotA;

xdot(5) = xdot(5)*((xdot(5)>0) || ~withinPad);

xdot = xdot';

